masterData <- read_csv("analysis_data.csv")

library(stargazer);library(sjPlot);library(gridExtra);library(ggrepel)
library(arm); library(scales);library(lme4)

set.seed(123)


##### 

# Extra regression part 1

# Table C1 Model with unanimous cases included


mControlsIdeologyLC_UNAN<-glmer(direction ~ scIdeolScore*unan+ 
                                  practiceSpecEconomic +practiceSpecCriminal+practiceSpecPublic+
                                  practiceSpecCivil+practiceSpecCommon+PriorExp+female+
                                  lcDispositionDirection +lcDissent+
                                  (1|ChiefJustice)+ (1|HCDBcaseId)+
                                  (1|primaryIssueSubArea)+(1|justice),
                                data=masterData, family = 'binomial',
                                control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))


tab_model(mControlsIdeologyLC_UNAN, show.est= T,show.re.var	=F,transform = NULL, show.r2	=F)




plot_model(mControlsIdeologyLC_UNAN, type = "int",terms='scIdeolScore', 
           colors = "Set2")+
  ggplot2::scale_y_continuous(limits = c(0,1))+ ggplot2::theme_classic()+
  ggplot2::scale_x_continuous(limits = c(0,1))+ggtitle("")+
  ggplot2::xlab("Liberal Ideology")+ggplot2::ylab('pr(Liberal vote)')
ggplot2::ggsave(filename = 'Figure_Unanimous_Interaction.png' ,height = 6, width = 8, device='png', dpi=600)

#showing that the null effect of unanimous cases does not depend on the area of law...

mControlsIdeologySpecPracticeVaryingSlopes <-glmer(direction ~  scIdeolScore*Area+
                                                     practiceSpecEconomic +practiceSpecCriminal+practiceSpecPublic+
                                                     practiceSpecCivil+practiceSpecCommon+PriorExp+female+
                                                     lcDispositionDirection +lcDissent+
                                                     (1|ChiefJustice)+ (1|HCDBcaseId)+
                                                     (1|primaryIssueSubArea)+(1|justice),
                                                   data=masterData[masterData$Area!="International and Maritime Law",], family = 'binomial',
                                                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
summary(mControlsIdeologySpecPracticeVaryingSlopes)


# interaction unanimous cases included (Table C2)
tab_model(mControlsIdeologySpecPracticeVaryingSlopes, show.est= T,show.re.var	=T,transform = NULL, show.r2	=F)





##### 



# Accounting for panel size C3



masterData$full_panel<-0
masterData$full_panel[masterData$numberJustices==7]<-1



mControlsIdeologyLC_ALL<-glmer(direction ~ scIdeolScore*full_panel+ 
                             practiceSpecEconomic +practiceSpecCriminal+practiceSpecPublic+
                             practiceSpecCivil+practiceSpecCommon+PriorExp+female+
                               lcDispositionDirection +lcDissent+
                             (1|ChiefJustice)+ (1|HCDBcaseId)+
                             (1|primaryIssueSubArea)+(1|justice),
                           data=masterData[masterData$unan==0,], family = 'binomial',
                           control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))


summary(mControlsIdeologyLC_ALL)


#Regression table


tab_model(mControlsIdeologyLC_ALL, show.est= T,show.re.var	=T,transform = NULL, show.r2	=F)

plot_model(mControlsIdeologyLC_ALL, type = "int", ci.lvl = NA ,terms='scIdeolScore', 
           colors = "Set2")+
  ggplot2::scale_y_continuous(limits = c(0,1))+ ggplot2::theme_classic()+
  ggplot2::scale_x_continuous(limits = c(0,1))+ggtitle("")+
  ggplot2::xlab("Liberal Ideology")+ggplot2::ylab('pr(Liberal vote)')
ggplot2::ggsave(filename = 'Figure_Full_Panel_Interaction.png' ,height = 6, width = 8, device='png', dpi=600)

##### 

# Table C4

## models replacing one type of pre appointment specialisation with another. 

mControlsIdeologySpecLC<-glmer(direction ~  scIdeolScore+ 
                                 LC_Spec_Economic +LC_Spec_Criminal+LC_Spec_Public+
                                 LC_Spec_Civil+LC_Spec_Common+LC_Spec_Procedure+PriorExp+female+
                                 lcDispositionDirection +lcDissent+
                           (1|ChiefJustice)+ (1|HCDBcaseId)+
                           (1|primaryIssueSubArea)+(1|justice),
                         data=masterData[masterData$unan==0,], family = 'binomial')

tab_model(mControlsIdeologySpecLC, show.est= T,show.re.var	=F,transform = NULL, show.r2	=F)



##### 

#Table C5

# Models with extra controls included

# running the 'Ideology' model several times with different sets of controls
# excluding all cases with fewer than 7 justices (full bench)



#ideology model

mControlsIdeology_1<-glmer(direction ~ scIdeolScore+ 
                                 practiceSpecEconomic +practiceSpecCriminal+practiceSpecPublic+
                                 practiceSpecCivil+practiceSpecCommon+PriorExp+female+
                                 lcDispositionDirection +lcDissent+
                                 (1|ChiefJustice)+ (1|HCDBcaseId)+
                                 (1|primaryIssueSubArea)+(1|justice),
                               data=masterData[masterData$unan==0,], family = 'binomial',
                               control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))


#ideology model with extra justice characteristics
mControlsIdeology_2<-glmer(direction ~ scIdeolScore+ 
                             practiceSpecEconomic +practiceSpecCriminal+practiceSpecPublic+
                             practiceSpecCivil+practiceSpecCommon+PriorExp+female+
                             lcDispositionDirection +lcDissent+
                             Experience+IsChiefJustice+
                             (1|ChiefJustice)+ (1|HCDBcaseId)+
                             (1|primaryIssueSubArea)+(1|justice),
                           data=masterData[masterData$unan==0,], family = 'binomial',
                           control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))




#ideology model with extra personal characteristics AND extra case facts

mControlsIdeology_3 <-glmer(direction ~ scIdeolScore+ 
                             practiceSpecEconomic +practiceSpecCriminal+practiceSpecPublic+
                             practiceSpecCivil+practiceSpecCommon+PriorExp+female+
                             lcDispositionDirection +lcDissent+
                             Experience+IsChiefJustice+
                             numberJustices+numFedGovParty+numStateGovParty+
                             Intervener+
                             (1|ChiefJustice)+ (1|HCDBcaseId)+
                             (1|primaryIssueSubArea)+(1|justice),
                           data=masterData[masterData$unan==0,], family = 'binomial',
                           control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))


#ideology model with extra personal characteristics AND extra case facts AND Experience interaction

mControlsIdeology_4 <-glmer(direction ~ scIdeolScore*Experience+ 
                              practiceSpecEconomic +practiceSpecCriminal+practiceSpecPublic+
                              practiceSpecCivil+practiceSpecCommon+PriorExp+female+
                              lcDispositionDirection +lcDissent+
                              IsChiefJustice+
                              numberJustices+numFedGovParty+numStateGovParty+
                              Intervener+
                              (1|ChiefJustice)+ (1|HCDBcaseId)+
                              (1|primaryIssueSubArea)+(1|justice),
                            data=masterData[masterData$unan==0,], family = 'binomial',
                            control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))

#Regression table
stargazer(mControlsIdeology_1,mControlsIdeology_2,mControlsIdeology_3,mControlsIdeology_4,
          out='Appendix Regression Table Ommitted Variable.html', align=T,style="ajps",
          star.cutoffs = c(.05, .01, .001)) 

#interaction plot
InteractionPlot<-plot_model(mControlsIdeology_4, type = "int", terms=c('Experience', 'scIdeolScore'), colors='bw')+
  ggplot2::ggtitle('')+
  ggplot2::scale_y_continuous(limits = c(0,1))+ ggplot2::theme_classic()+
  ggplot2::scale_x_continuous(limits = c(0,1))+
  ggplot2::xlab("Liberal Ideology")+ggplot2::ylab('pr(Liberal vote)')

InteractionPlot

ggplot2::ggsave(InteractionPlot,filename = 'supplementary_Regression_Interaction_Figure.png' ,height = 5, width = 7, device='png', dpi=600)

##### 


# Table C6

# running identical models to those in the main manuscript, 

# But replacing the chief justice reign RE with the natural court RE C6 


mIdeology_NC<-glmer(direction ~    scIdeolScore +
                   (1|ChiefJustice)+ (1|HCDBcaseId)+
                   (1|naturalCourtDecision)+(1|justice),
                 data=masterData[masterData$unan==0,], family = 'binomial',
                 control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))

mControlsIdeology_NC<-glmer(direction ~    scIdeolScore +
                           practiceSpecEconomic +practiceSpecCriminal+practiceSpecPublic+
                           practiceSpecCivil+practiceSpecCommon+PriorExp+female+
                           (1|naturalCourtDecision)+ (1|HCDBcaseId)+
                           (1|primaryIssueSubArea)+(1|justice),
                         data=masterData[masterData$unan==0,], family = 'binomial',
                         control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))

mControlsIdeologyLC_NC<-glmer(direction ~ scIdeolScore+ 
                             practiceSpecEconomic +practiceSpecCriminal+practiceSpecPublic+
                             practiceSpecCivil+practiceSpecCommon+PriorExp+female+
                               lcDispositionDirection +lcDissent+
                             (1|naturalCourtDecision)+ (1|HCDBcaseId)+
                             (1|primaryIssueSubArea)+(1|justice),
                           data=masterData[masterData$unan==0,], family = 'binomial',
                           control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))


summary(mControlsIdeology_NC)
summary(mControlsIdeologyLC_NC)

#Regression table
stargazer(mIdeology_NC, mControlsIdeology_NC,mControlsIdeologyLC_NC,
          out='Appendix Regression Table Natural Courts.html',align=T,style="ajps",
          star.cutoffs = c(.05, .01, .001))




###Table C7

source("Table 6.R")

## model 5 bayesian model using new data and bootstrapped standard error

dlist<- list(ideol = comparison$meanBoot, 
             ideol_SD = comparison$sdBoot, 
             meanDirection = model_sum$meanDirection, 
             N = nrow(model_sum), 
             numberCases = model_sum$numberCases, 
             numberSuccesses = model_sum$numberSuccesses,
             
             #controls
             female = model_sum$female, 
             PriorExp = model_sum$PriorExp, 
             econ_exp = model_sum$practiceSpecEconomic, 
             crim_exp = model_sum$practiceSpecCriminal, 
             publ_exp = model_sum$practiceSpecPublic, 
             civi_exp = model_sum$practiceSpecCivil, 
             comm_exp = model_sum$practiceSpecCommon, 
             
             
             lcDispositionDirection = model_sum$lcDispositionDirection, 
             lcDissent = model_sum$lcDissent
)




model_5 <- ulam(
  alist(
    
    #regression specification
    #numberSuccesses ~ dbinom(numberCases, p),
    
    meanDirection~dnorm(mu, sigma),
    mu <- a +bI*ideol_true[i] + 
      b_econ_exp*econ_exp + b_crim_exp*crim_exp + b_publ_exp*publ_exp + b_civi_exp*civi_exp + b_comm_exp*comm_exp + 
      bPriorExp*PriorExp + bFem*female +
      b_lcDispositionDirection*lcDispositionDirection + b_lcDissent*lcDissent, 
    
    #accounting for uncertainty in estimated ideology
    ideol ~ dnorm(ideol_true, ideol_SD),
    vector[N]:ideol_true ~dnorm(0,1),
    
    #standard priors
    bI ~dnorm(0,0.5),
    b_econ_exp ~dnorm(0,0.5),
    b_crim_exp ~dnorm(0,0.5),
    b_publ_exp ~dnorm(0,0.5),
    b_civi_exp ~dnorm(0,0.5),
    b_comm_exp ~dnorm(0,0.5),
    bPriorExp ~dnorm(0,0.5),
    bFem ~dnorm(0,0.5),
    b_lcDispositionDirection ~dnorm(0,0.5), 
    b_lcDissent ~dnorm(0,0.5),
    
    a~dnorm(0,0.2),
    
    sigma~dexp(1)
  ),data = dlist, chains=4, cores=4, iter= 2000)


# Model 5 coumns of table 7 
model_5_table<-precis(model_5,prob = 0.95)[,1:4]

